The following is a demo of a turning point analysis (Jones, 2017) to determine the change from enrichment to dilution in 4 Iowa catchments. Each catchment has daily N concentration and water discharge data from 2015 through 2019 collected by Iowa State University.
Map developed by Dr. Peiyu Cao.
The raw data is stored in an Excel file, with each catchment separated on a different sheet. This function imports the data for each site and calculates the area-normalized discharge (water yield). This is then stored in a list for future access.
# import data
site_fc <- function(n){
# one catchment per sheet
c1 <- read_excel("Data/Peiyu/Catchments_raw.xlsx", sheet = n)
# rename columns
c1 <- select(c1,
date = Date,
Q_m3_day = `Inflow (m3/day)`,
c_mg_L = `[NO3] in (mg/L)`,
Interpolation) %>%
filter(Interpolation == "F") %>% # only keep raw data
mutate(site = n,
name = excel_sheets("Data/Peiyu/Catchments_raw.xlsx")[n],
year = year(date),
month = month(date),
Q_norm_mm = case_when(
name == "KS" ~ (Q_m3_day/3067200)*1000, # divide Q by drainage area in m2
name == "RS" ~ (Q_m3_day/4806900)*1000, # and multiply by 1000
name == "LP" ~ (Q_m3_day/2596500)*1000, # to convert to mm
name == "WW" ~ (Q_m3_day/2188800)*1000
)) %>%
select(-Interpolation)
return(c1)
}
site_list <- lapply(1:4, site_fc)Once the data is imported, we made a custom function that separates the data by year. It then calculates the Parametric Quantile Regression (PQR) using the Ricker function as described in Jones, 2017.
The output of the PQR model returns the coefficients (a, b, c) for the equation. These coefficients can then be used to calculate the turning point for each year. The custom function is then applied to the 4 catchments.
Note: discharge is filtered to retain values above zero (in order to perform log functions). Higher discharge values are not yet filtered.
PQR_fc <- function(n){
year <- c(2015:2019)
name <- excel_sheets("Data/Peiyu/Catchments_raw.xlsx")[n]
# make fc to filter data by year
year_fc <- function(y){
sub <- as.data.frame(site_list[[n]]) %>%
filter(year == y, Q_norm_mm > 0) %>% # log won't work on negative values
select(date,
x = Q_norm_mm,
y = c_mg_L)
return(sub)
}
# separate all years
year_list <- lapply(year, year_fc)
# make a new fc to apply PQR model for each year
model_fc <- function(i){
df <- as.data.frame(year_list[[i]])
# plug equation into model
model_rq <- rq(y ~ log(x) + x, data = df)
coef <- model_rq$coefficients
return(coef)
}
# apply to all years, individually
model_list <- lapply(1:length(year), model_fc)
# calculate turning point based on model coefficients
tp_fc <- function(i){
a <- model_list[[i]][[1]]
b <- model_list[[i]][[2]]
c <- model_list[[i]][[3]]
# Q in mm, area-normalized
tp <- -b/c
return(tp)
}
tp_list <- lapply(1:length(year), tp_fc)
# the output will be one table with the tp values (by year)
tp_vector <- c(tp_list[[1]],tp_list[[2]],tp_list[[3]],tp_list[[4]],tp_list[[5]])
df <- as.data.frame(cbind(year, name))
df <- df %>% cbind(tp_vector)
colnames(df) <- c("year", "name", "tp_Q_mm_norm")
return(df)
}
# all turning points by catchment-year data point
PQR_list <- lapply(1:4, PQR_fc)
PQR <- rbindlist(PQR_list)Turning points (tp) as area-normalized discharge:
When the custom PQR function is run with the data for the 4 catchments, some of the turning points have negative values. This can be further explored by visualizing the modeled curves.
Building off from the custom PQR function, a new one was made to automate the generation of the following plots.